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ABSTRACT 

Three  methods  for  generating  outcomes  on  multivariate 
normal  random  vectors  are  presented.   A  comparison  is  made 
to  determine  which  method  requires  the  least  computer 
execution  time  and  memory  space  when  utilizing  the 
IBM  36O/67.   All  methods  use  as  a  basis  a  standard  Gaussian 
random  number  generator.   Results  of  the  comparison  study 
indicate  that  the  method  based  on  triangular  factorization 
of  the  covariance  matrix  generally  requires  less  memory 
space   and  computer  time  than  the  other  two  methods. 
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I.   INTRODUCTION 

Frequently  it  is  desired  to  generate  a  sample  of  vectors 
derived  from  a  specified  multivariate  normal  distribution. 
For  example,  this  need  arises  in  simulation  studies  in  which 
the  model  contains  stochastic  variates  that  are  distributed 
according  to  a  multivariate  normal.   Since  a  computer  is 
frequently  used  to  produce  the  sample  vectors,  the  memory 
space  and  computer  time  requirements  become  important  cost 
factors  when  considering  various  methods  for  "generating" 
random  vec'  )rs . 

The  problem  of  "generating"  random  vectors  was  discussed 
in  19^8  by  Herman  Wold  [Ref.  1] .   Wold  describes  a  method 
for  construction  of  samples  from  a  multivariate  normal 
distribution.   Wold  used  a  method  based  on  triangularization 
of  the  covariance  matrix.   In  1962,  Ernest  M.  Scheuer  and 
David  S.  Stoller  wrote  a  paper  on  the  generation  of  multi- 
variate normal  random  vectors  [Ref.  2] .   They  considered  a 
method  based  on  matrix  equations  and  a  method  based  on  con- 
ditional distributions.   The  first  method  uses  a  technique 
similar  to  that  of  Wold's.   Their  conclusion  was  that  the 
mathematics  involved  in  the  first  method  was  simpler  than 
that  of  the  second  method  so  the  generation  of  random 
vectors  is  best  accomplished  by  the  matrix  equation 
technique . 

The  objective  of  this  paper  is  to  examine  the  above  two 
methods,  hereafter  called  "triangularization"  and 


"conditioning"  respectively,  and  to  discuss  a  third  method 
which  is  referred  to  as  the  "rotation"  method.   Computer 
programs  for  each  method  are  presented  and  a  comparison  of 
memory  space  requirements  and  execution  times  on  the 
IBM  36O/67  computer  is  made  for  the  three  methods.   The 
format  for  making  this  comparison  study  is  similar  to  that 
used  by  M .  E.  Muller  in  his  paper  on  univariate  normal 
generators  [Ref .  31  • 

The  following  notation  is  used  in  this  study:   X  denotes 
a  random  variable  whereas  x  is  a  real  variable.   A  vector 
is  denoted  by  a  bar  under  the  letter.   Listed  below  are  two 
examples  of  the  above  notation: 


X  = 


X 


1 


X. 


X 
L  PJ 


x  = 


x 
L  PJ 


In  the  discussion  of  the  three  methods,  generators  for 
N(0_,£)  distributions  are  considered,  where  £  is  a  general 
(positive  definite  symmetric)  matrix.   There  is  no  loss  of 
generality  in  assuming  the  mean  vector  to  be  the  zero 
vector  since  a  random  vector  Y  distributed  N(u,E)  may  be 
generated  by  first  generating  X  distributed  N(0_,Z)  and  then 
taking  Y  =  X  +  u. 


II.   METHODS 

A.   ROTATION  METHOD 

Given  a  positive  definite  symmetric  p  x  p  matrix  E,  the 
function 

f(x;E)  =  (27r)-^p  |  E  |  ~h   e~^'Z   -  (1) 

is  a  p-variate  normal  density.   A  random  vector  X  having 
this  distribution  has  mean  0_  and  variance  covariance  matrix 
E  [Ref.  k] .   Denote  the  distribution  law  corresponding  to 
the  density  function  (1)  by  N(p_,E).   If  X  (with  p  components) 
is  distributed  N(0,E),  then  Y  =  CX  is  distributed  M(0,CEC), 
if  C  is  nonsingular  [Ref.  H] .   The  development  of  the 
rotation  method  is  based  on  this  fact.   If  a  p  x  p  matrix  C 
is  found  such  that  CEC  yields  a  diagonal  matrix,  then 
elements  on  the  diagonal  of  CEC'  are  the  variances  of  the 
corresponding  components  of  Y.   Since  the  covariance  values 
Cpff  diagonal  elements  of  CEC')  are  zero,  the  components  of 
Y_  are  independent. 

For  a  symmetric  matrix  E,  there  exists  an  orthogonal 
matrix  C  such  that  CEC  =  D  where  D  contains  the  eigen 
values  of  E  on  its  diagonal  [Ref.  5] •   Associated  with  the 
eigen  values  are  eigen  vectors  which  make  up  the  rows  of 
the  matrix  C.   Since  Y  =  CX  is  N(0,CEC)>  and  there  exists 
a  C  such  that  CEC  is  the  diagonal  matrix  D,  it  follows 
that  the  components  of  Y  are  independent  normal  variates 
with  zero  mean  and  variance  values  as  given  by  the 


corresponding  eigen  values  of  Z.   These  variates  may  be 
generated  using  a  univariate  normal  random  number  generator 
A  method  of  generating  the  random  vector  Y  is  thus  easily 
established.   To  obtain  a  generator  for  the  random  vector 
X,  a  transformation  is  made  using  X  =  C~  Y.   Since  C  is 
an  orthogonal  matrix,  this  simplifies  to  the  equation 
X  =  C'Y.   Standard  matrix  multiplication  of  the  orthogonal 
matrix  C'  and  sample  vectors  generated  on  the  random 
vector  Y  yields  the  desired  random  samples  of  the  vector  X. 

B.   CONDITIONING  METHOD 

A  multivariate  probability  density  function  f  may  be 
written  as: 


f(x13...,xn)  =  f(x1|x2,. . . ,xn)  f(x2|x3,. . . ,xn) 

. . .  f(x   - |x  )  f(x  )  (2) 

n-11  n     n 


where  f(x.|x....)  denotes  a  conditional  density  function. 
J  '  i  _ 

If  X   can  be  generated  from  f(x  ),  then  X   ,  may  be  gener- 
n        &  n  '       n-1   J  to 

ated  by  conditioning  on  X  .   X  _?  may  then  be  generated  by 

conditioning  on  X   and  X   n .   In  a  similar  manner,  all  the 
to     n      n-1  3 

components  of  X  may  be  generated.   In  what  follows,  this 
general  approach  is  applied  to  multivariate  normal 
distributions . 

Given  that  the  vector  X  is  partitioned  into  two  sub- 
vectors  X    and  X    and  that  the  covariance  matrix  of  X 
is  correspondingly  partioned  into  £-1-.,  ^ip»  an<^  ^??}  then 
it  can  be  shown  [Ref  *J]  that  the  conditional  distribution 
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.,  v(l)   .     v(2)     (2)  .        ,   ...       _    r  -1   (2) 
of  X    given  X    =  x     is  normal  with  mean  E.„  E„?   x 

and  covariance  matrix  E,,  -  E.„  E  ~   E„  .   Several  defini- 
tions are  introduced  at  this  time  as  they  are  required  in 
what  follows. 
Let 


(k)  _ 


k,k 


k+l,k 


k,k+l 


k+l,k+l 


k,n 


k+l,n 


0      i 
n,k 


n,k+l 


and  let  E . .   denote  the  cofactor  of  a.,  in  I    :  similarly 
E.  .  denotes  the  cofactor  of  a.  .  in  E.   Consider  now  the  mean 
and  variance  of  the  distribution  determined  by 
f  (x,  |  x,  ,..,...  ax  ).   Suppose  the  covariance  matrix  E  is 
partioned  as  follows : 


E  = 


Zll=  [ok,k] 


21 


k+1 

,k 

k+2 

,k 

n,k 

12    k,k+l     k,k+2        k,nJ 


22 


°k+l,k+l  °k+l,k+2  "•  ak+l,n 


ak+2,k+l   ak+2?k+2  ""  ak+2,n 


n,k+l 


n,k+2 


...  a. 


(3) 
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The  variance  of  f  (x,  |  x,  ,  ,  .  .  .  ,x  )  is  given  by: 
V(Xklxk+l'"  ,,xn)  =  Zll  ~  Z12  E22   E21 


(h) 


To  simplify  this  expression,  the  following  theorem  is  used 

[Ref.  H]  :       |E|  =  ]£■>-,    -    £-,p  ^22   ^p-i  I-'  I  ^pp  I  •   The  square 
matrix  Z„„  is  nonsingular  since  the  covariance  matrix  Z  is 
positive  definite.   The  variance  (4)  can  be  written  as: 


00 


V(Xk'xk+ls"  *,xn)  "  |v(k+i) 


(5) 


The  mean  of  f (xk| xk+1 , . . . ,xn)  is  given  by 


E(Xklxk+l'- ' ,,Xn)  "  Z12  Z22   -' 


Direct  substitution  of  components  of  the  partitioned  covari- 
ance matrix  (3)  into  the  latter  expression  yields: 


E(X,  |x.  ........ x  )  =  [a.    ,  ,  n  a,  .  ,  „  ...  a,  ](Z  ) 

k'  k+ls    '  n     k,k+l   k,k+2      k,n 


Lk+1 


lk+2 


n 


A  simplification  of  this  expression  can  be  made  as  follows 
First,  examine  (E      )   .   By  using  standard  matrix 
algebra,  it  can  be  shown  that: 
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(z(k+l)} 


-1 


(k+1) 


.(k+1) 
'k+1, k+1 


-E 


(k+1) 
k+l,k+2 


-E 


(k+1) 
k+2,k+l 


.(k+1) 

k+2,k+2 


,_1sk+l+nE(k+l)       ,_1xk+2+nz(k+l) 
k+l,n  k+2,n 


(-if 


+k+ly(k+l) 
n,k+l 


(_1)n+k+2E(k+l) 
n,k+2 


(      )2n  j.  (k+1) 

n,n 


Next,  perform  the  matrix  multiplication 


(ak,k+l   ak,k+2   •'•   ak,n]  [Z 


(k+1) 


-1 
]   . 


This   yields    a   vector  of   dimension    1   x    (n   -   k  +    1).      The 

transpose   of  this    vector  is    the   product    of  the    factor 
1 


Tk+T7 


and 


"[  +  E(k+1)      a 

1      k+1, k+1  k, k+1 


[-E(k+1)      a 

1      k+2,k+l   k,k+l 


E(k+1)      a 
k+l,k+2   k,k+2 


E(k+1)      a 
k+2,k+2   k,k+2 


(_     k+l+nz(k+l)  k+2+nz(k+l) 

n,k+l    k,k+l  n,k+2 


k+l,n   k,n      1 


( _1^+2+nI(k+l)  „ 

k+2,n  k,n      d. 


(_1)2ns(k+l)        , 

n,n       k,n        3 


(6) 
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Multiply  the  vector  (6)  by  [x,  -.,..., x  ]  and  the  following 
results  are  obtained: 


E(Xk|xk+1,...,xn)  =~n^Ty7  .[(al)xk+l+(a2)xk+2  ••'+  (a3)xn] 

where  (a.)  indicates  the  first  element  of  the  vector  (6), 
(a~)  indicates  the  second  element,  etc.   Upon  close  inspec- 
tion, it  is  apparent  that  (a,)  is  equivalent  to  -E^.n  u.- 

In  the  same  manner,  (a~)  is  -Z, .  L   ,  .  .  ..,  and  (a-,)  is 

'    2       k+2,k'     '       3 

-£   ; .   The  mean  of  f(x,  |x,  .,...., x  )  can  thus  be  represented 
n,k  k1  k+1'    '  n  ^ 


as 


E(XklXk+r*  •*jXn)  ""  iz(k+l)i  [Ek+l,kxk+l  +  Ek+2^xk+2--*+Zn,kXn] 

(7) 
The  use  of  these  expressions  in  generating  the  components 
of  X  is  now  described.   Consider  the.  transformation 
X,  =  o.    Y.  +  y,  ,  k  =  1,2,...  n,  where  the  Y's  are  indepen- 
dent normal  variates  with  zero  mean  and  unit  variance. 
Using  expressions  (5)  and  (7)  for  o      and  y,  respectively, 
this  transformation  may  be  given  explicitly  as: 


X  =  y   /a  (8) 

n    n   n,n 

_  r  u(k)i  1%      ,  i     rr(k)        (k)  y       +y^y~\ 

Ak   [|E(k+l)|   *k   ij-Ck+l)  I  [_k+l,kAk+l  Zjk+2,kAk+2*  *  *   n,k  nj 


(9) 
for  k  =  l,...,n-l.   Use  of  the  equations  (8)  and  (9)  there- 
fore provides  a  method  of  generating  samples  of  a  random 
vector  X,  using  univariate  generators   for  the  Y's. 
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C.   TRIANGULARIZATION  METHOD 

The  covariance  matrix,  designated  as  Z,  is  a  symmetric 
positive  definite  matrix.   It  can  be  written  as  Z  =  TT' 
where  T  is  a  lower  triangular  matrix  [Ref .  6] , 


all   a12 


a21   a22 


nl   n2 


In 


2n 


nn 


— I 

^™ 

hi 

0 

0 

hi 

fc21 

.  .  .   t   , 

nl 

t21 

U   r-y  *-\                   •      i 

0 

0 

t22 

...   tn2 

nl 

t     0          .  . 

n2 

.       t 
nn 

0 

0 

• 
• 
• 

*  *  *  fcnn 

X  is  obtained  by  making  the  non-singular  transformation 
X  =  TY  where  Y  is  distributed  N(0_,I).   By  previous  remarks, 
X  is  distributed  H(0_,TT').   The  matrix  T  can  be  obtained 
from  Z  by  using  the  so  called  "square  root  method"  [Ref.  6] 
Equations  for  t. .  in  terms  of  a. .  are  as  follows: 


t.n  =  c.n/ 


il  "  uilVa 


11 


il 


V     i"1 


=1   IP 


;   1  <  i  <  n 


1  <  i  <  n 


1        J"1 

; .  .  =  [a..-Z   t.   t .  ]  :   l<j<i<n 

ij   ffjj    id   p=i   iP  dP   ' 


t.  .  =  0 
id 


;  i  <  d   n 


The  vector  Y  consists  of  independent  normal  components  with 
zero  mean  and  unit  variance.   These  may  be  generated  using 
a  univariate  normal  random  number  generator.   Standard 
matrix  multiplication  of  T  and  samples  of  Y  yields  samples 
of  the  random  vector  X. 
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III.   METHOD  EVALUATION 

Computer  runs  were  obtained  using  each  of  the  three 
methods,  with  various  size  covariance  matrices.   Data  con- 
cerning these  runs  is  compiled  in  Tables  I,  II,  and  III. 
A  brief  description  of  the  procedure  employed  to  obtain 
execution  time  and  memory  space  requirements  for  the  tables 
follows.   The  computer  programs  can  be  subdivided  into  three 
basic  portions.   The  first  portion  consists  of  the  data 
that  must  be  entered  into  the  computer.   The  second  portion 
consists  of  all  required  steps,  prior  to  the  use  of  the 
basic  univariate  random  number  generator,  necessary  to 
implement  one  of  the  three  methods  described  in  this  study. 
The  third  portion  consists  of  the  standard  Gaussian  random 
number  generator  and  those  steps  necessary  for  generating 
samples  of  the  random  vector,  X,  such  that  a  matrix  is 
filled  with  one  thousand  random  vectors.   Two  times  were 
considered  when  evaluating  the  programs:  the  "set-up"  time 
and  the  "repetition"  time.   These  times  correspond  to 
execution  time  for  the  second  and  the  third  portion  respec- 
tively of  the  computer  program.   Similarly,  the  memory 
space  requirements  shown  in  the  tables  make  use  of  two 
numbers.   The  first  number  provides  an  indication  of  the 
amount  of  space  required  for  the  computer  programs  for 
each  of  the  three  methods.   This  number  includes  the  sub- 
routine EIGEN  in  the  rotation  method  and  the  function  GRN 
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in  all  cases.   The  second  number  Is  the  total  space 
required  for  the  program  plus  any  external  function  used 
such  as  square  roots,  absolute  values,  exponentials,  input 
and  output  devices,  etc.   The  second  number  varies  as  to 
the  computer  system  being  used  and  in  this  case  the  computer 
was  the  IBM  36O/67.   The  headings  on  two  columns  of  Tables  I, 
II,  and  III  are  labeled  "Specific  Conditioning"  and 
"Specific  Triangularization"  respectively.   These  refer  to 
computer  programs  written  specifically  for  2x2  and  3x3 
covariance  matrices,  and  are  adaptations  of  the  general 
programs . 

Based  on  the  data  contained  in  Tables  I,  II,  and  III, 
the  best  method  to  use  for  generating  random  vectors  from 
covariance  matrices  dimensioned  greater  than  3x3  appears 
to  be  the  triangularization  method.   This  method  requires 
less  memory  space  and  considerably  less  execution  time  than 
the  other  two  methods.   If  the  user  is  interested  in 
generating  samples  from  a  bivariate  or  trivariate  normal 
distribution,  then  either  a  specific  conditioning  or  triang- 
ularization method  can  be  used.   Both  methods  require 
about  the  same  amount  of  execution  time  and  their  space 
requirements  are  essentially  the  same. 

All  methods  use,  as  a  basis,  the  univariate  normal 
random  number  generator.   Therefore,  in  order  to  establish 
a  lower  bound  on  times  required  in  generating  normal 
random  vectors,  the  time  required  to  generate  normal  random 
numbers  was  obtained.   These  numbers  should  serve  as  lower 
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bounds  on  the  times  required  to  generate  a  sample  of  one 
thousand  repetitions  of  a  10,  5,  3,  and  2  element  vector 
respectively.   The  lower  bounds  obtained  in  this  way  are: 

10,000  numbers  -  1. 067*10  seconds 
5,000  numbers  -  0.53209  seconds 
3,000  numbers  -  0.32302  seconds 
2,000  numbers  -  0.20332  seconds 

A  specialized  covariance  matrix,  (identity  matrix), 
was  used  as  an  input  to  all  three  methods  to  provide  the 
least  cumbersome,  mathematically  speaking,  of  all  matrices 
that  would  be  encountered.   The  results  of  this  test  are 
tabulated  in  Table  III.   Each  method  maintained  its  overall 
ranking  with  respect  to  memory  space  and  time  requirements. 
The  rotation  method  displayed  a  sharp  decrease  in  set-up 
time  which  is  reasonable  as  eigen  values  are  easier  (faster) 
to  compute  in  this  case.   In  general,  each  of  the  three 
methods  required  essentially  the  same  amount  of  time  for 
1000  repetitions  using  a  typical  covariance  matrix  as  for 
the  specialized  covariance  matrix. 
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TABLE  I 
MEMORY  SPACE  REQUIREMENTS 

The  first  number  is  the  amount  of  core  space   required 
for  the  computer  program  and  the  second  number  is  the 
total  amount  of  core  space  required  in  execution  of  the 
program.   This  second  number  includes,  for  example, 
external  functions  and  input-output  devices.   Each 
space  requirement  is  in  bytes. 


Matrix 
Size 

Conditioning 

Triangular- 
ization 

Rotation 

2x2 

11,792/37,168 

9,528/33,864 

11,952/36,288 

3  x  3 

15,880/41,256 

13,568/37,904 

16,008/40,344 

5x5 

24,160/^9,536 

21,712/46,048 

24,192/48,528 

10  x  10 

45,424/70,800 

42,336/66,672 

44,992/69,328 

-  -  -  - 

Specific  Cond. 

Specific  Triang. 

-  -  -  - 

2x2 

8,752/33,088 

8,744/33,080 

-  -  -  - 

3  x  3 

12,976/37,312 

12,960/37,296 

-  -  -  - 
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TABLE  II 
EXECUTION  TIME  (TYPICAL  COVARIANCE  MATRIX) 

The  first  number  is  the  set-up  time,  and  the  second 
number  is  the  repetition  time.   Each  execution  time  is 
in  seconds. 


Matrix 
Size 

Conditioning 

Triangular- 
ization 

Rotation 

2x2 

.00273/. 33985 

.00106/. 33800 

.00496/. 37666 

3x3 

.00741/. 55383 

.00145/. 52411 

.01398/. 62023 

5x5 

.03975/1.06335 

.00247/. 99554 

.06117/1.23357 

10  x  10 

.70306/2.75560 

.00795/2.52280 

.41652/3.72289 

_.  _  _  _ 

Specific  Cond. 

Specific  Triang. 

-.  -  -  - 

2x2 

.00091/. 25246 

.00101/. 27487 

-  -  -  - 

3x3 

.00101/. 39530 

.00109/. 39572 

-  -  -  - 
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TABLE  III 
EXECUTION  TIME  (SPECIALIZED  COVARIANCE  MATRIX  I) 

The  first  number  is  the  set-up  time,  and  the  second 
number  is  the  repetition  time.   Each  execution  time  is 
in  seconds. 


Matrix 
Size 

Conditioning 

Triangular- 
ization 

Rotation 

2x2 

.00278/. 3^299 

.00122/. 33509 

.00171/. 36376 

3x3 

.00793/. 55869 

.00161/. 5^223 

.00234/. 60250 

5x5 

.0321^/1.06025 

.00260/. 98449 

.00403/1.19545 

10  x  10 

.33795/2.69337 

.00733/2.42645 

.01162/3.41616 

-  _  «  _ 

Specific  Cond. 

Specific  Triang. 

*-  -  -  - 

2x2 

.00104/. 25215 

.00078/. 24921 

_  _  _  _ 

3x3 

.00135/. 37903 

.00096/. 38795 

-  -  -  - 
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APPENDIX  A 
DISCUSSION  OP  THE  COMPUTER  PROGRAMS 

All  three  computer  programs  were  written  in  Fortran  IV. 
They  are  not  optimal,  i.e.,  improvements  to  the  programs 
are  possible.   The  results  of  these  improvements  could 
change  the  value  of  memory  space  requirements  and  program 
execution  times,  but  should  not  change  the  comparative 
ranking  among  the  three  methods. 

Two  univariate  normal  random  number  generators  that 
were  available  for  use  in  the  present  study  were  GAUSS  and 
GRN.   The  subroutine  GAUSS  is  part  of  the  IBM  scientific 
subroutine  package  and  is  based  on  the  mean  value  theorem. 
The  function  GRN  was  compiled  by  the  Naval  Postgraduate 
School  staff,  and  is  based  on  the  G.  Marsaglia  technique 
[Ref .  7] .   The  execution  time  of  GRN  is  approximately  ten 
times  faster  than  GAUSS,  and  for  this  reason  it  was  used 
in  all  the  programs.   The  rotation  method  uses  the  sub- 
routine EIGEN  which  is  also  available  as  part  of  the  IBM 
scientific  subroutine  package.   Three  subroutines,  (DTERM, 
REDMAT,  and  COFACT),  were  developed  for  use  in  the  con- 
ditioning method.   The  subroutine  DTERM  computes  the  deter- 
minant of  an  n-square  matrix.   The  subroutine  COFACT 
removes  the  elements  from  the  i   row  and  j^n  column  of  an 
n-square  matrix  (Z),  and  provides  the  (n  -  l)-square  matrix 
that  is  required  in  the  computation  of  the  cof actor  Z. .. 
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The  subroutine  REDMAT  reduces  a  matrix  to  the  desired  size. 

(k) 
It  is  capable  of  providing  the  user  with  a  Z     from  an 

n  x  n  dimension  matrix  where  k  <_   n. 

Each  method  was  evaluated  using  four  different  size 
matrices  as  listed  in  Tables  I,  II,  and  III.   These  were 
selected  as  being  typical  of  the  range  of  matrix  sizes  that 
might  be  encountered  in  practice.   However,  the  programs 
were  written  in  general  terms  so  there  are  no  program 
imposed  restrictions  on  matrix  size. 

As  mentioned  previously,  specific  computer  programs 
were  written  for  the  2x2  and  the  3x3  covariance  matrices 
using  the  conditioning  and  triangularization  method.   This 
was  done  in  an  attempt  to  evaluate  the  possible  reduction 
in  memory  space  and  execution  time  under  those  obtained 
using  the  general  programs  with  input  dimensions  of  2  and 
3<   Of  the  three  methods  considered,  these  two  appeared  to 
be  adaptable  to  simpler,  more  concise  programs  for  these 
small  dimensions.   Writing  specific  programs  for  covariance 
matrices  dimensioned  greater  than  3x3  would  become  dif- 
ficult, and  a  general  program  would  probably  have  to  be 
used.   No  attempt  was  made  to  write  a  specific  program  for 
the  rotation  method  since  computation  of  eigen  values  and 
vectors  are  involved,  and  there  were  no  apparent  means  of 
reducing  the  computation  time  under  that  of  using  the  general 
subroutine  EIGEN. 
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The  computer  programs,  used  in  generating  multivariate 
normal  random  vectors  for  each  of  the  three  methods,  are 
included  as  a  portion  of  this  study. 


2h 


COMPUTER  PROGRAMS 


THIS  PROGRAM  UTILTZFS  THE  ROTATION  MPTHOn  TQ  GENERATE 
MULTIVARIATE  NORMAL  RANDOM  VECTORS. 


DIMFNSION  7(3,3),  R(3,3),  XA(3,1000),  R  ( f> )  ,  Y(3),  TO) 
READ  (5,1)  N 

1  FORMAT  (15) 
DO  2  I  =  1,N 

2  PtAO  ( c,3)  (7(  I , J)  ,  J  -  1 ,N) 

3  FORMAT  (3FP.0) 
MV  =  0 

K  =  1 

DO  5  J  =  1,N 
DO  4  I  *  It J 
B  (  K  )  =  Z  (  I  ,  J  ) 
^     K  =  K  +  1 

5  CONTINUE 
LB  =  K  -  1 

CALL  EIGEN  (B,R,N,MV) 

MA  =  1 

JA  =  1 

DO  6  LA  =  1,N 

Y(  JA)  =  B(MA) 

Y(  JA)  =  SQRT(  Y(  JM  ) 

JA  =  JA  +  1 

6  MA  =  MA  +  JA 
IDUMMY  =  0 

DO  10  MT  =  1,1000 

DO  7  JP  =  1,N 

X  =  (GRN( IDUMMY) )*Y( JP) 

7  T(  JP )  =  X 

DO  9  I  -    1VN 
SUM  =  CO 
DO  8  J  =  1,N 

8  SUM  =  SUM  4-  R(I  t  J)  *  T(  J) 

9  XA(  I,MT)  =  SUM 

10  CONTINUE 
STOP 

END 
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THIS  PROGRAM  UTILIZE  THE  TP I ANGUl A P T Z AT  I PN  MFTHPD 
TO  GENERATE  MULTIVARIATE  NORMAL  RANDOM  VECTPRS. 


7 
8 


9 
10 

11 

12 

13 
14 


DIME 
READ 
FORM 
READ 
FORM 


I  =  1  ,  N) 


DO 
C( 
IF 
JE 
DP 
C( 


L 

If 
( 

If 

CONT 
CO  NT 
DO  1 
IF  ( 
JT  = 
DO  7 
SUM 
JD  = 
DO  6 
SUN 
CONT 
C(  I, 
CONT 
SUM 
JF  = 
DO  9 
SUM 
CONT 
C(  I, 
CONT 
IDUM 
DO  1 
DO  1 
X  = 
T(  JC 
DO  1 
SUMM 

DO  1 
SUMM 
XV(  J 

CONT 
CONT 
STOP 
END 


NSICN  Z(3,3),  C(3f3)f  XVQOOO,?), 

(  ^,1)  N 
AT  (IP) 

(  5f 2)  ( (Z(  I , J)  ,  J  =  1 t N)  , 

4T  (3FP.0) 

I  =  ltN 
1)  =  7(1 ,1) /SORT  (7(1,1)) 

I  +  1  .GT.  M)  GO  TO  & 

I  +  1 

JR  =  J',N 
J  P  )  =  0.0 
INLF 
INUF 

0  I  =  2fN 
T  -  3)  2,5,5 

I  -  1 
J  =  2,JT 
=  0.0 
J  -  1 
ka  =  i,jn 

=  SUM  +  C( I  ,KA)  *  C( J,KA) 

INUE 

J  )  =  (  Z  (  I  ,  J  )  - 

INUE 

=  0.0 

I  -  1 

K  =  ltJP 
=  SIM  +  C( I  ,K) 
INUE 

1  )  =  SORT  (7(1,1) 

INUE 

MY  =  0 

**    JP  =  1,1000 

1  JC  =  1,N 
GRN( I  DUMMY) 
)  =  X 
3  IT  =  1,N 

=  0.0 

2  J  =  1,TT 
=  SUMM  +  C( IT, J) 

PfIT)  =  SUMM 

INUE 

INUE 


(3) 


SUN) /C( J, J) 

*  C  (  I  ,  K  ) 
SUM) 


*  T(  J) 
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THIS    PROGPAM    UTILI7CS    THC    CONDITIONAL    MFTHOD    TP 
GENERATE    MULTIVARIATE"    NORMAL    RANDOM    VICTORS. 


DIMENSION    Z(3,3)t     M3,3),    B(3f?)f    Y(3),     R(3),    T(3) 
1XA( 1000,3) ,    ZJ( 3,3) 
READ    («j.l)     N 

1  FORMAT     (IA) 

READ    (5t2)     UZ(I,J),     J    =    l.N'li    I     =    ltN) 

2  FORMAT     (3F8.0) 
M    =    N    -    1 

DO    6    L    =    1,M 

K    =    N    -    L 

LT    =    N    -    K    +    1 

CALL    REDMAT    (Z,N,LT,K,A) 

IR    =    K    +     1 

DO    3    IV    =    IR,    N 

NV    =    LT    -     1 

IF    =     IV    -    K    +    1 

IG    =     1 

CALL    COFACT    ( LT  ,  A ,8 , NV , I F , I G) 

CALL    DTFPM    (NV,    B,    D,     NV) 

IT    =     IF    +    IG 

3  ZJ(K  ,  I  V)    =    ( (-1)**IT)     *    0 
CALL    DTFPM    (LT,A,D,LT) 
DNUM    =    D 

DO    A    I S    =    1,N 
DO    4    JS    =    1,N 
A(  IS, JS)     =    0.0 
B( ISf JS  )    =    0.0 

4  CONTINUE 

CALL    REDMAT    ( Z ,  f !  ,N V , I R  ,  A ) 

CALL    DTEPM(NV,A,D,NV) 

T(K)     =    D 

DO    5    IC    =    1,N 

DO    5    JC    =    1,N 

A(  ICJC  )    =    0.0 

5  CONTINUE 

6  R(K)     =    SORT(DNUM/T(K)  ) 
AD    =    SORT(Z(N,N) ) 
IDUMMY    =    0 

DO    10    NOP    =    1,1000 
DO    7    NO    =    l,N 
X    =    GRN( IDUMMY) 

7  Y(NO)    =    X 
XA(NOP,N)     =    AD    *    Y(N) 
LA    =    N    -    1 

DO    9    NORM    =    1,LA 

SUMN    =    0.0 

K    =    N    -    NOPM 

LEE    =    K    +    1 

=    LFF,N 

K,NSYL)     *    XA(NOP,NSYL)    +    SUMN 

=    R(  K)     *    Y(K)     -     ( SUMN/T(K) ) 


DO    8    NSYL 

8 

SUMN    =    ZJ( 

9 

XA(NOP  ,K) 

10 

CONTINUE 

STOP 

END 
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SUBROUTINE    DTFRM    (N,A,P,M) 

PIMFNSION    A(M,M) 

OD    =    1.0 

00    1°    L    =    lfN 

KP    =    0 

7    =    0.0 

00     12    K    =    L,N 

IF     (Z    -    ABS( A(K ,L) ) )     11,12,12 

11  Z    =    ABS(A( K,L) ) 
KP     =    K 

12  CONTINUF 

IF    (L    -    KP)     13,15,15 

13  DO     1&    J    =    L,N 
Z    =    A(L,J) 

A  (  L  ,  J  )    =    A(KPtJ) 

14  A(KP,J)    =    Z 
DD    =    -    DO 

15  IF    (L    -    N)     16,20,20 

16  LP1    =    L    +    1 

DO    19    K    =    LPl    ,    N 

IF    ( A( K,L)  )     17,19,17 

17  RATIO    =    A(K,L) /A(L,L) 
00    lft    J    =    LP1,M 

IB  A(K,J)     =    A(K,J)    -    RATIO    *    A(L,J) 

19  CONTINUF 

20  00    21K=1,N 

21  OD    =    DD    *    A( K,K) 
D    =    DD 

RETURN 
END 
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2  a 

25 


SUB 

DIM 

I    = 

DO 

IF 

I    = 

J    = 

DO 

IF 

J    = 

A(  I 

CON 

CON 

RFT 

END 


ROUT 
ENS  I 

0 
23    L 
(L     . 

!    + 

0 
22    K 
<K     . 

J    + 
,J> 
TINU 
TINU 
URN 


INE    REDMAT    ( Z,N,LZ,IB,A) 
ON    Z(N,N) ,    A(L7,LZ) 

=    1,N 
LT.     IB)    GO    TO    23 
1 

=    1,N 
LT.    IB)    GO    TO    22 

1 
=    Z(L,K) 


SUBROUTINE    COFACT     ( NA , Z ,B , NV, I K , IL ) 

DIMENSION    Z(NA,NA),    B(NV,NV) 

I    =    0 

nn    7^    k    =    i . ma 


K    = 
.EG. 
+    1 


DO    25 

IF     (K 

I    =    I 

J    =    0 

DO    21 

IF    (L 

J    =    J 

B(  I  ,J)     = 

CONTINUE 

CONTINUE 

RETURN 

END 


1,NA 
IK) 


GO    TO    2  5 


L    =    1,NA 
.EQ.    IL) 
1 
Z(KtL) 


GO    TO    24 


28 


BIBLIOGRAPHY 


1.  Wold,  H. ,  Tracts  for  Computers,  No.  XXV,  Cambridge 

Press,  1955. 

2.  Scheuer,  E.  M.  and  Stoller,  D.  S.,  "On  the  Generation 

of  Normal  Random  Vectors,"  Technometrics ,  v.  4, 
no,  2,  p.  278-281,  May  1962. 

3.  Muller,  M.  E. ,  "A  Comparison  of  Methods  for  Generating 

Normal  Deviates  on  Digital  Computers,"  Journal  of 
the  Association  for  Computing  Machinery"^  v~.    6 ,  no.  3 » 
p.  376-3B3,  July  1959. 

4.  Anderson,  T.  W. ,  An  Introduction  to  Multivariate 

Statistical  Analysis,  p.  17-19,  App .  A,  Wiley,  i960. 

5.  Graybill,  F.  A.,  An  Introduction  to  Linear  Statistical 

Models,  p.  5,  v.  I,  McGraw-Hill,  196I. 

6.  Graybill,  F.  A.,  Introduction  to  Matrices  with  Applica- 

tions in  Statistics,  p.  298,  Wadsworth,  1969 . 

7.  Marsaglia,  G.,  "A  Fast  Procedure  for  Generating  Normal 

Random  Variables,"  Communications  of  the  ACM,  v.  7, 
no,  1,  p.  4-10,  January  1964. 

8.  Howe,  J.  E. ,  The  Generation  of  Random  Numbers  from 

Various  Probability  Distributions,  Masters  Thesis, 
Naval  Postgraduate  School,  Monterey,  1965. 


29 


INITIAL  DISTRIBUTION  LIST 

No.  Copies 

1.  Defense  Documentation  Center  2 
Cameron  Station 

Alexandria,  Virginia   2231^ 

2.  Library,  Code  0212  2 
Naval  Postgraduate  School 

Monterey,  California   939^0 

3.  Department  of  Operations  Analysis,  Code  55       1 
Naval  Postgraduate  School 

Monterey,  California   939^0 

h.      Associate  Professor  D.  R.  Barr,  Code  55  bn       1 
Department  of  Operations  Analysis 
Naval  Postgraduate  School 
Monterey,  California   939^0 

5.   LCDR  N.  L.  Slezak,  USN  1 

Box  3^ 
Milligan,  Nebraska   68^06 


30 


UNCLASSIFIED 


Security  Classification 


DOCUMENT  CONTROL  DATA  -R&D 

{Security  classification  of  title,    body  of  abstract  and  indexing  annotation  must  be  entered  when   the  overall  report  Is   classified) 


1.   originating    ACTIVITY   (Corporate  author) 

Naval  Postgraduate  School 
Monterey,  California   939^0 


2a.   REPORT    SECURITY    CLASSIFICATION 

Unclassified 


Zb.    GROUP 


3      REPOR  T     TITLE 


A  COMPARISON  OF  METHODS  FOR  GENERATING 
MULTIVARIATE  NORMAL  RANDOM  VECTORS 


4.   DESCRIPTIVE   NOTES  (Type  of  report  and,  in  cf  us  ive  dates) 

Master's  Thesis;  September  1Q70 


5.   AUTHORISI  (First  name,  middle  initial,  last  name) 


Norman  Lee  Slezak,  Lieutenant  Commander,  United  States  Navy 


6   REPOR  T  D A  TE 


September  1970 


7a.     TOTAL    NO.    OF    PAGES 


30 


76.    NO.    OF    REFS 


ta.    CONTRACT    OR    GRANT    NO. 


b.    PROJECT    NO. 


9a.    ORIGINATOR'S    REPORT    NUMBER(S) 


So.   OTHER   REPORT   NO(5>  (Any  other  numbers    that  may  be  assigned 
this  report) 


10.    DISTRIBUTION    STATEMENT 


This  document  has  been  approved  for  public  release  and  sale; 
its  distribution  is  unlimited. 


II.    SUPPLEMENTARY    NOTES 


12.    SPONSORING   MILITARY    ACTIVITY 


Naval  Postgraduate  School 
Monterey,  California  939^0 


13.  ABSTRAC  T 


Three  methods  for  generating  outc 
normal  random  vectors  are  presented, 
to  determine  which  method  requires  th 
execution  time  and  memory  space  when 
IBM  360/67.   All  methods  use  as  a  bas 
random  number  generator.   Results  of 
indicate  that  the  method  based  on  tri 
of  the  covariance  matrix  generally  re 
space  and  computer  time  than  the  othe 


omes  on  multivariate 

A  comparison  is  made 
e  least  computer 
utilizing  the 
is  a  standard  Gaussian 
the  comparison  study 
angular  factorization 
quires  less  memory 
r  two  methods , 


DD,F, 


0RM  1473 

nov  es  I  *t  I  sj 

S/N    0101 -807-681 1 


(PAGE     1) 


31 


UNCLASSIFIED 


Security  Classification 


A-91408 


UNCLASSIFIED 


Security  Classification 


key    wo  RDS 


MULTIVARIATE   NORMAL    DISTRIBUTION 
RANDOM   NUMBER   GENERATOR 
NORMAL   DISTRIBUTION 
MULTIVARIATE   NORMAL   GENERATOR 


ROLE  W  T 


ROLE  W  T 


P  .TvM..1473  'BACK 

0101  -807-6821 


ROLE  WT 


32 


UNCLASSIFIED 


Security  Classification 


A-3t  409 


225  15 


Thesis 

S5709 

C.l 


L20094 

Slezak 

A  comparison  ot 
methods   for   generatmc 


mu 


iltivariate  normal 
.  random  vecto^s2  515 


Thesi  s 

120094 

S5709 

Slezak 

-  1 

A  comparison  of 

methods  for  generating 

multivariate  normal 

random  vectors. 

3  2  '68  001  006  17  4 

•   KNOX  I 


